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ABSTRACT 

We explore the effect of magnetorotational turbulence on the dynamics and concentrations of boul- 
ders in local box simulations of a sub-Keplerian protoplanetary disc. The solids are treated as particles 
each with an independent space coordinate and velocity. We find that the turbulence has two effects 
on the solids. 1) Meter and decameter bodies are strongly concentrated, locally up to a factor 100 
times the average dust density, whereas decimeter bodies only experience a moderate density increase. 
The concentrations are located in large scale radial gas density enhancements that arise from a com- 
bination of turbulence and shear. 2) For meter-sized boulders, the concentrations cause the average 
radial drift speed to be reduced by 40%. We find that the densest clumps of solids are gravitationally 
unstable under physically reasonable values for the gas column density and for the dust-to-gas ratio 
due to sedimentation. We speculate that planctcsimals can form in a dust layer that is not in itself 
dense enough to undergo gravitational fragmentation, and that fragmentation happens in turbulent 
density fluctuations in this sublayer. 

Subject headings: instabilities — MHD — planetary systems: formation — planetary systems: proto- 
planetary disks — turbulence 



1. INTRODUCTION 

Planets are believed to form from micrometer-sized 
dust grains that grow by collisio nal sticking in p rotoplan- 
etary ga s discs (iSafronoy lllQfi.Ql see reviews bv iLissauerl 
1993 and lBeckwith et all200(fl . Once the bodies reach a 
size of around one kilometer, the growth to Moon-sized 
protoplanets and later real p lanets is achieved by gravi- 
tationally induced collisions ijThommes et alJ|2008l) Al- 
though significant progress has been made in the un- 
derstanding of the initial conditions of grain growth 
ijHenning et alJl2005fl . we nevertheless do not yet have 
a complete picture of how the solids grow 27 orders of 
magnitude in mass to form kilometer-sized planetesimals. 

Growth by coagulation can take place when there is a 
relative speed between the solids. Various physical effects 
induce relative speeds at different grain size scales. This 
allows for a definition of distinct steps in the growth from 
micrometer dust grains to meter-sized boulders in a tur- 
bulent protoplanetary disc. Microscopic dust grains gain 
their relative speed due to Brownian motion. This pro- 
cess forms relatively comp act cluster-cluster aggregates 
ijDominik fc Tielenslll997li . The speed of the Brownian 
motion falls rapidly with increasing grain mass, and so 
the time-scale for building up larger compact bodies this 
way becomes prohibitively large, compared to the life- 
time of a protoplanetary disc. 

When Brownian motion is no longer important, the 
relative speed is dominated by the differential vertical 
settling in the disc. The vertical component of the cen- 
tral star's gravity causes the gas to be stratified. Dust 
grains do not feel the pressure gradient of the gas and 
thus continue to fall towards the mid-plane with a veloc- 
ity given by the balance between vertical gravity and 
the drag force. Larger grains fall faster than smaller 
grains due to the size-dependent coupling to the gas (ac- 
tually bodies that are so massive that they are starting 
to decouple from the gas will rather move on inclined or- 
bits relative to the disc, i.e. perform damped oscillations 



around the mid-plane). As they fall, they are thus able to 
sweep up smaller grains in a process that is qualitatively 
similar to rainfall in the Earth's atmosphere. Upon ar- 
rival at the mid-plan e, the largest s olids can reach sizes 
of a few centimeters ( Safronov 1969). These bodies have 
grown as compact particle-cluster aggregates with a high 
porosity. 

Turbulent gas motions cause the sedimented solids 
to diffuse away from the mid-plane llCuzzi et a,lj ri993: 
iDubrulle et"aTTll995() . where they can meet and collide 
with a reservoir of microscopic grains. These tiny grains 
still hover above the mid-plane because their sedimen- 
tation time-scale is so long that turbulent diffusion can 
keep them well-mixed with the gas over a large verti- 
cal extent. Turbulence also plays a role for equal-sized 
macroscopic bodies by inducing a relative collision speed 
that is much larger than the Brownian moti on contribu- 
tion ijVolk et al.lll980HWeidenschillinelll984|) . 

When estimating the outcome of an interaction be- 
tween macroscopic bodies, the issues of collision physics 
must be taken into account. For relative speeds above a 
certain threshold, the bodies are lik ely to break up when 
they collide rather t han to stick (jChokshi et alJ 119931: 
Blu m fc WurTi]l2000]) . This is a problem for macroscopic 
bodies where the sticking threshold is a few meters per 
second. Fragmentation caused by high-speed encounters 
continuously replenish the reservoir of microscopic dust 
grains. These can then be swept up by the boulders that 
are lucky enough to avoid critical encounters. However, 
the sweeping up of smaller dust grains by a macroscopic 
body has its limitations w hen the relative spe ed exceeds 
some 10 meter per second ijWurm et al.l l2001). At larger 
relative velocities of up to a hundred meters per second, 
which are likely to occur due to the high speed of larger 
bodies, the small grains will erode the boulder. 

The time evolution of the size-distribution of 
solids can be calculated by solvi ng the coagu- 
lation equation numerically (e.g. iWetheriil 119901 
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Weidenschillindll997t iSuttner fe Yorkell2001() . Recently, 
Dullcmond & Dominik ( 2005) performed numerical sim- 
ulations of the coagulation for realistic disc environ- 
ments. Starting with micrometer-sized grains only, they 
find that a narrow peak of 0.1-10 meter-sized boulders 
can form in 10 4 -10 5 years, when fragmentation is ig- 
nored. On the other hand, in a more realistic situa- 
tion high speed impacts lea d to fragmentation. Here 
IDullemond fc Dominikl lj2005|) find that once the size dis- 
tribution reaches the meter regime, still around 75% of 
the mass is maintained in microscopic bodies, which are 
the fragments of larger bodies that have been destroyed 
in collisions. This picture is given some credit by the fact 
that microscopic dust grains are observed in protoplan- 
etary discs of millions of years of age, whereas the time- 
scale for depleting grains of those sizes is only around 
1,000 years in the absence of fragmentation. 

Besides the problem of getting macroscopic bodies 
to stick, meter-sized boulders quickly drift radially in- 
ward toward the central star due to their aerodynamic 
friction with the gas in a typical sub-Keplerian disc 
ijWeide nschillingl Il977j) . The drift time-scale can be as 
short as 100 years. To avoid evaporation in the inner 
disc or in the central star, the bodies must grow by least 
an order of magnitude in size (three orders of magnitude 
in mass) in a time shorter than this! 

A possibility to overcome the growth obst acles was 
suggested independently by ISafronovl l)1969|) and by 
Goldr eich fc Wardl l)1973(L The general idea is that boul- 
ders sediment towards the mid-plane and form a particle 
sublayer that undergoes a gravitational instability, form- 
ing the planetesimals in a spontaneous event (gelation) 
rather than by continuous growth (coagulation). The 
weakest point in this model is that it requires a laminar 
disc in order to work. Even a tiny amount of turbulence 
in the disc will prevent the boulders from an efficient 
sedimentation towards the mid-plane, and the insta - 
bility will never occur l|Weidenschilling fc Cuz^ll [19931. 
Thus disc turbulence had always to be avoided in order 
to allow for self-gravity assisted planetesimal formation. 
However, even in a completely laminar disc, the settled 
dust induces a vertical shear in the gas rota tion profile 
l|Weidenschilhndll980t iNakagawa et alJH986f) . This can 
be unstable to a Kelvin-Helmholtz instability. The sub- 
sequent Kelvin-Helmholtz turbulence puffs up the dust 
layer so that the densities needed for a gravitational in- 
stability are usually not achieved, unless a dust-to-gas 
ratio ma ny times higher than the solar composition is 
adopted iYoudin fc ShiJl20T)l . 

Nevertheless, solids can reach sizes of around one 
meter without the help of self-gravity In this size 
regime the gradual decoupling from the gas motion en- 
ables the bodies to move independently from the gas. 
This can cause them to be trapped in turbulent fea- 
tures of the gas flow. An important theoretical dis- 
covery is that meter-sized bo ulders are concentrated in 
gaseous anticyclonic vortices |Barge fc Sommerial 119951 
IChavanisI 120001: IJohansen et al]l2004j) . Inside such vor- 
tices the dust density can locally be enhanced to val- 
ues sufficient either for enhanced coagulation or even 
for gravitational fragmentation. Also the radial drift 
of part icles trapped in the vortices is si gnificantly re- 
duced ijde la Fuente Marcos fc Bargdl200l|) . Theoreti- 
cal attention has furthermore been given to the trap- 



ping of dust grains in high pressure regions. Since dust 
grains do not feel pressure forces, any pressure-supported 
gas structure must cause dust grai ns to move i n the 
direction of the pres s ure gr adient lIKlahr fc LirJ l200ll 
iHaghighipour fc Boss! 1200% iKlahr fc LirJ l2005lK Re- 
cently, |2i££_£L^D (|2QQ4[) demonstrated that this can lead 
to large concentrations (a density increase of up to a fac- 
tor 50) of meter-sized boulders in the high density spi- 
ral arms of self-gravitating discs. The same mechanism 
can drain millimeter-sized dust grains from the under- 
dense regions around a protoplanet that is not massive 
enou gh to open a gap in the gaseou s component of the 
disc i|PaardekooDer fc Mellemall2004lK 

Giant long-lived vortices may form in pro- 
toplanetary disc due to a baroclinic instability 
ijKlahr fc Bodenheimerl I2003D . but the conditions 
for the barocli nic instability in protoplanetary discs are 
still not clear (Klahr 2004). Magnetorotational turbu- 
lence (MRI) on the other hand is expected to occur in 
all discs where t h e ionization fraction is sufficiently high 
llGammiei 119961 iFromang et al.l l2002t iSemenov et alJ 
|2004|) . A search for dust conce ntrations in magnetorota- 
tional turbulence was done by iHodgson fc Brandenburg! 
(1998) who found n o apparent concent rations. On the 
other hand, recently I Johansen fc Klahrl (2005, hereafter 
referred to as JK05) found evidence for centimeter-sized 
dust grains being trapped in short-lived turbulent eddies 
present in magnetorotational turbulence. That work 
was, however, limited by the fluid description of dust 
grains, i.e. the friction time must be much shorter than 
the orbital period, and could not handle grains larger 
than a few centimeters. 

In this paper, we expand the work done in JK05 by 
putting meter-sized dust particles, represented by real 
particles rather than by a fluid, into magnetorotational 
turbulence. We show tha t magnetorotational turbulence 
(Balbus & Hawlcy 1991) is not actually an obstacle to 
the self gravity-aided formation of planetesimals, but 
rather can be a vital agent to produce locally gravita- 
tional unstable regions in the solid component of the disc 
when the average density in solids would not allow for 
fragmentation. This process is very similar to the gravo- 
turbulent frag mentation of molecular clouds into pro- 
tostellar cores (|Klessen et al.ll200fl iPadoan fc Nordhmdl 
2004). 

2. DYNAMICAL EQUATIONS 

For the purpose of treating meter-sized dust boulders 
we ha ve adapted the Pencil Code 1 (see also lBrandenburel 
2003) to include the treatment of solid bodies as particles 
with a freely evolving (2, y, z)-coordinate on top of the 
grid. This is necessary because the mean free path of 
the boulders, with respect to collisions with the gas 
molecules, is comparable to the scale height of the disc. 
Thus the dust component can no longer be treated as a 
fluid, but must be treated as particles each with a freely 
evolving spatial coordinate Xi and velocity vector Vi . In 
other words, it is no longer possible to define a unique 
velocity field at a given point in space for the particles, 
because they keep a memory of their previous motion. 
Friction only erases this memory for small grains. 

1 The code is available at 
http : //www.nordita. dk/sof tware/pencil-code/ 
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TABLE 1 
Simulation parameters 



Run 


JV 


Lq; X Ly X Lz 


n x X n y X ra z 


n 


Qqt { 


p 


At 


A 


2 x 10 6 


1.32 x 1.32 x 1.32 


64 x 64 x 64 


7.6 


1.0 


-0.04 


100 


B 


2 x 10 6 


1.32 x 1.32 x 1.32 


64 x 64 x 64 


7.6 


0.1 


-0.04 


100 


C 


2 x 10 6 


1.32 x 1.32 x 1.32 


64 x 64 x 64 


7.6 


10.0 


-0.04 


100 


D 


2 x 10 6 


1.32 x 1.32 x 1.32 


64 x 64 x 64 


7.6 


1.0 


0.00 


100 


E 


2 x 10 6 


1.32 x 5.28 x 1.32 


64 x 256 x 64 


1.9 


1.0 


-0.04 


24 


F 


2 x 10 6 


1.32 x 10.56 x 1.32 


64 x 512 x 64 


1.0 


1.0 


-0.04 


16 



Note. — First column: name of run; second column: number of particles; third 
column: size of the box measured in scale heights; fourth column: grid dimension; 
fifth column: number of particles per grid cell; sixth column: friction time; seventh 
column: global pressure gradient parameter; eighth column: number of orbits that 
the simulation has run. 



2.1. Drag Force 

The particles are coupled to the gas motion by a drag 
force that is proportional to the velocity difference be- 
tween the particles and the gas, 



/drag = - — (vi-u) . 

Tf 



(1) 



Here u is the gas velocity at the location of particle i and 
rf is the friction time. The friction time depends on the 
solid radius a, and the solid density p, as 



Tf = 



a 2 , p. 



min(a.c s , \v)p ' 



(2) 



where v is the molecular viscosity of the gas, c s is the 
sound speed and p is the gas density. This expression 
is valid when the particle speed is m uch lower than the 
sound speed \ Weidenschillind 1 1 977|) . Using the kinetic 
theory expression for viscosity v — c s A/2, where A is the 
mean free path of the gas molecules, the friction time can 
be divided into two regimes: the Epstein regime is valid 
when a, < 9/4A. Here the mean free path of the gas 
molecules is longer than the size of the dust grain, so the 
gas can not form any flow structure around the object. 
The friction time is proportional to the solid radius in 
this regime. In the Stokes regime, where a. > 9/4A, 
a flow field forms around the object. Now the friction 
time is proportional to solid radius squared, so the object 
decouples faster from the gas with increasing size. For 
an isothermal and unstratified disc, one can treat the 
friction time Tf as a constant. The distinction between 
the Epstein and the Stokes regime is then only important 
for translating the friction time into a solid radius (see 
end of this section) . 

To determine the gas velocity in equation at the po- 
sitions of the particles, we use a three-dimensional first- 
order interpolation scheme, using the eight grid corner 
points surrounding a given particle. For multiprocessor 
runs the particles can move freely between the spatial 
intervals assigned to each processor using MPI (Message 
Passing Interface) communication. 

2.2. Disc Model 



We consider a protoplanetary disc in the shearing sheet 
approximation, but for a disc with a radial pressure gra- 
dient d\nP/d\nr = a (or P oc r a ). In the shear- 
ing sheet approximation this gradient produces a con- 
stant additional force that points radially outwards (be- 
cause the pressure falls outwards). Making the vari- 
able transformation In p — > In p + (l/ro)ax, the stan- 
dard isothermal shearing sh eet equation of motion (e.g. 
Gol dreich fc Tremainelll978h gets an extra term, 



du 
~dt 



+ (u-V)u = -2f2 xu + 3f2„x~c 2 V\np- 



'0 



(3) 

The terms on the right-hand-side of equation © are the 
Coriolis force, the centrifugal force plus the radial gravity 
expanded to first order, and the two terms representing 
local and global pressure gradient. The coordinate vector 
(x, y,z) is measured from the comoving radial position 
ro from the central source of gravity, with x pointing 
radially outwards and y along the Keplerian flow. At 
r = ro the Keplerian frequency is Qq. The shearing 
sheet approximation is valid when all distances are much 
shorter than ro. The balance between pressure gradient, 
centrifugal force and gravity is given for a sub-Keplerian 
rotation of the disc, 



(4) 



where the first term on the right-hand-side is the purely 
Keplerian rotation profile, while the second (constant) 
term is the adjustment due to the global pressure gradi- 
ent. We now measure all velocities relative to the sub- 
Keplerian flow using the variable transformation u — > 
u + Uq. This changes equation @ into 



9u , 

w + («.v)«. 



•u(°>^ = /(«)-c?Vlnp. (5) 



Here the last term on the left-hand-side represents the 
advection due to the rotation of the disc relative to the 
center of the box (which moves on a purely Keplerian 
orbit). The function / is defined as 



/(«) 



2f2ou y 1 

— \QqU x 

. , 



(0) 



When making the same variable transformation in the 
equation of motion of the dust particles, there is however 
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TABLE 2 

Results 



Run 


O Tf 


p 


max(n) 




(lam) 


(7 






t7z 


A 


1.0 


-0.04 


81.3 


— 0.0123 


-0.020 


0.0222 


0.0162 


0.0105 


0.0077 


B 


0.1 


-0.04 


32.6 


— 0.0034 


-0.004 


0.0139 


0.0064 


0.0101 


0.0052 


C 


10.0 


-0.04 


77.5 


-0.0042 


-0.004 


0.0170 


0.0115 


0.0094 


0.0062 


D 


1.0 


0.00 


56.5 


-0.0003 


0.000 


0.0225 


0.0165 


0.0106 


0.0078 


F 


1.0 


-0.04 


50.3 


-0.0132 


-0.020 


0.0204 


0.0149 


0.0093 


0.0066 


E 


1.0 


-0.04 


50.3 


-0.0132 


-0.020 


0.0194 


0.0140 


0.0086 


0.0061 



Note. — First column: name of run; second column: friction time; third column: global 
pressure gradient parameter; fourth column: maximum particle density in units of the 
average density; fifth column: radial velocity averaged over space and time; sixth column: 
predicted radial drift in a non-turbulent disc; seventh to tenth columns: velocity dispersion 
averaged over space and time. Averages are taken from 5 orbits and beyond. Grid cells 
with or 1 particles have been excluded for the calculations of velocity dispersions. 



no global pressure gradient term to balance the extra 
Coriolis force imposed by the sub-Keplerian part of the 
motion, so the result is 



dvj 

at 



i 



/(«,•) (Vi -u) + 



Tf 



2 1 - 

Cg — ax . 

ro 



(7) 



The modified Coriolis force / appears again because of 
the presence of Xi(t) in ito. The last term on the right- 
hand-side reflects the head wind that the dust feels when 
it moves through the slightly sub-Keplerian gas. The 
reason that the term appears in the radial component of 
the equation of motion is that all velocities are measured 
relative to the rotational velocity of the gas. A dust 
particle moving at zero velocity with respect to the gas 
thus experiences an acceleration in the radial direction. 

The explicit presence of ro in equation Q is non- 
standard in the shearing sheet. It may seem that the 
term vanishes for ro — > oo. But this is actually not the 
case, since the natural timescale of the disc, f^Q 1 , also 
depends on ro, so that at large radii there is an immense 
amount of time at hand to let the tiny global pressure 
gradient force work. One can quantify this statement by 
dividing and multiplying by the scale-height H in the 
last term of equation (JJJ to obtain the result 



dv. 



Of 



c s (2 — a . 
ro 



(8) 



Here H/r = £ is the ratio of the scale height to the 
orbital radius, a quantity that is below unity for thin 
discs. Depending on the temperature profile of a disc, 
the typical value of £ is between 0.001 and 0.1. We define 
the pressure gradient parameter f3 as (3 = a£. 

For the simulations, we adopt the following dynamical 
equations for gas velocity u, magnetic vector potential 
A, gas density p, particle velocities i>j and particle coor- 
dinates Xf. 

^- + (u-V)u + u y °^ = /(«) - c s 2 V lnp 



1 



J xB + /„(«,/*) (9) 



dA (m OA „ 



- f2 A y x + fJA) (10) 



dy 
dvj 
dt 



■f(vi) + c s n /3x 

1 



-(Vi - u) 



dxj 
~~dt 



(11) 

(12) 
(13) 



The functions /„, /„ and /b are hyperdiffusivity terms 
present to stabilize the finite difference numerical scheme 
of the Pencil Code. This is explained in more detail in 
JK05. We shall ignore the effect of the global pressure 
gradient on the dynamics of the gas, since for ^ <C 1 the 
increase in density due to the global gradient is much 
smaller than the average density in the box. Thus we 
set simply Uy = —3/2f2ox. We also ignore the con- 
tribution from the global density on the Lorentz force 
term in equation and the advection of global density 
in equation (|ll|l . Furthermore we do not include verti- 
cal gravity in the simulations. This means that we solve 
exactly the same equations for the gas as in JK05, i.e. 
without radial pressure stratification. The radial drift 
of solids then originates exclusively from the dynamical 
equations of the particles. 

We solve the dynamical equations (|9")l- (|13[) for vari- 
ous values of the friction time and of the box size. The 
typical resolution is 64 3 for a box size of 1.32H on all 
sides. A similar setup was used in JK05 to calculate 
the turbulent diffusion coefficient of dust grains in mag- 
netorotational turbulence. In the present work we ex- 
pand the model by letting 2,000,000 particles represent 
the dust grains. Thus the dust component is typically 
represented by approximately 8 particles per grid cell. 
We set the strength of the radial pressure gradient by 
the parameter f3 — —0.04. This would represent e.g. a 
disc with a global pressure gradient given by a = — 1 and 
a scale-height-to-radius ratio o f £ = 0.04, which is typ- 
ical f or a solar nebula model ( Wcidenschilli n~fc Cuzzj 
1993). We consider friction times of f2 T { = 0.1,1,10. 
The translation from friction time into grain size depends 
on whether the friction force is in the Epstein or in the 
Stokes regime, but the two drag laws yield quite similar 
grain sizes in the transition regime. Thus, at the radial 
location of Jupiter in a typical protoplanetary disc, the 
friction time corresponds to grains of approximately 0.1, 
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Fig. 1. — The number of particles in the densest grid cell as a function of time for run A (meter-sized boulders). The maximum density 
is generally around 20 times the average, but peaks at above 80 times the average particle density. The insert shows a magnification of the 
time between 50 and 51 orbits. 




Fig. 2. — The number of particles in the densest grid cell as a function of time, here for runs B (decimeter-sized boulders) and C 
(decameter-sized boulders). The first shows only very moderate overdensities, whereas the latter is similar in magnitude to run A (meter- 
sized boulders), but with broader peaks. 



1 and 10 meters in size. 

The simulation parameters are given in Table ^ We 
let the boulders have random initial positions from the 
beginning and let them start with zero velocity. 



3. PARTICLE CONCENTRATIONS 

In Fig. ^ we plot the number of particles in the dens- 
est grid cell as a function of time for run A (meter-sized 
boulders, see Table QJ. The average number of particles 



per grid cell is 7.6. Evidently there is more than 100 
particles in the densest grid cell at most of the times, 
and at some times the number is even above 600. This is 
more than 80 times the average dust number density. In 
Fig. |21 we plot the maximum particle density for runs with 
r2o r f = 0.1 (run B, gray curve) and J?o r f = 10 ( run C, 
black curve) . The decimeter-sized boulders are obviously 
not as strongly concentrated as the meter-sized boulders, 
whereas the decameter-sized boulders have concentra- 
tions that are similar in magnitude to run A. The mea- 
sured values of the maximum particle density for all the 



Johansen et al. 



0.97 



I/In 



1.03 0.0 



I d /(e Io) 





-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

x/(c s ^) 



-0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 

x/(c s Qo 1 ) 



Fig. 3. — Gas column density E (left panel) and dust column density J7<j (right panel). The gas column density only varies by a few 
percent over the box, but still a slightly overdense region is seen near the center of the box. The dust column density in the same region is 
up to 5 times the average dust column density. 



runs can be found in Table [5] 

To examine whether some structures in the gas density 
arc the source of the high particle densities, we plot in 
Fig-dlthe column densities of gas £ and of dust particles 
Zd at a time of 50.9 orbits for run A. The gas column 
density varies only by a few percent over the box, since 
the turbulence is highly subsonic, but a region of moder- 
ate overdensity is seen around the middle of the box. The 
dust column density is very high in about the same re- 
gion as the gas overdensity, around a factor of five higher 
than the average dust column density in the box, so dust 
particles have moved from the regions that are now un- 
derdense into the overdensity structure near the center 
of the box. 

We explore the radial density structure of the gas and 
the dust in the box in more detail in Fig. 0] Here the 
azimuthally averaged gas and dust column densities are 
shown as a function of radial position x and time t mea- 
sured in orbits. Apparently large scale gas density fluctu- 
ations live for a few orbits at a constant radial position 
before decaying and reappearing at another radial po- 
sition. The fluctuation strength is less than 1% of the 
average density. The dust density shows strong peaks 
at the locations of the gas density maxima. The ex- 
planation for this correlation is as follows. Locations 
of maximal gas density are also local pres sure maxima. 
Such pressure maxima can trap du st grains ijKlahr fc Linl 
l200lt iHaehighipour fc Boss1l2003|) as they are locations 
of Keplerian gas motion. The inner edge of a pressure 
maximum must move faster than the Keplerian speed be- 
cause the pressure gradient mimics an additional radial 
gravity. At the outer edge of a radial pressure enhance- 
ment the outwards-directed pressure gradient mimics a 
decreased gravity, and the gas must move slower than 
the Keplerian speed. Dust grains do not feel the pressure 
gradient and are thus forced to move into the pressure 
bump. In our simulations the radial gas overdensities 



have a typical lifetime at a given radial position on the 
order of a few orbits. When the gas overdensity even- 
tually disappears, the particle overdensity is only slowly 
getting dissolved, and the particles drift and concentrate 
towards the location of the next gas overdensity. The gas 
density structure in the azimuthal and vertical directions 
does not show a similar density increase, and as expected 
there is also no significant concentration of particles with 
respect to these two directions. The density fluctuations 
thus have the form of two-dimensional sheets. 

In Fig. [3] we plot the maximum density experienced by 
a 200 particle subset of the 2,000,000 particles during the 
100 orbits. The distribution function £(n) is defined as 
the fraction of particles that have been the center of a 
number density of at least n over the size of a grid cell. 
The curves clearly show how large the concentrations are. 
For decimeter-sized boulders, 95% of them have experi- 
enced a 5 times increase in dust density, whereas only 
around 2% have been part of a 10 times increase. For 
meter-sized particles, 70% have been part of a 10 times 
increase in dust density, and 1% even took part in a 20 
times increase. Particles of decameter-size had more than 
10% taking part in a 30 times increase of d ust density. 
This i s very similar to the concentrations that lRice et Id] 
( 2004) find in the spiral arms of self-gravitating discs. 

In Fig. correlations between gas flow and particle 
density are shown for run D (without global pressure 
gradient). Here we have taken data at every full or- 
bit, starting at 5 orbits when the turbulence has satu- 
rated, and calculated the average particle density in bins 
of various gas parameters. We also plot the spread in 
the particle density in each bin. The top two panels 
show the correlation with two components of the vortic- 
ity u) = V x u. There is some correlation between ver- 
tical vorticity component and the particle density, but 
the spread in each vorticity bin is larger than the av- 
erage value. The correlation indicates that some trap- 
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Fig. 4. — Azimuthally averaged gas and dust column densities as a function of radial position relative to the center of the box x and time 
t. Black contour lines are shown at gas density fluctuations of 0.5% from the average value. Large scale density fluctuations are seen to 
have lifetimes on the order of a few orbits before moving to other radial positions. The dust column density peaks strongly at the locations 
of the maximal gas column density. 

particle concentrations stay together even after the gas 
feature which created them has decayed or moved to an- 
other location. The limited life-time of the concentrating 
features weakens the measured correlation with the gas 
flow. 

The lower two panels of Fig. show the correlation 
with divergence of pressure gradient flux and with gas 
density. In a steady flow, particles accelerate towards 
an equilibrium velocity where the drag force is in bal- 
ance with the other forces working on the particles. The 
equilibrium velocity is 

v = np- 1 (VP- J x B) = T f F. (14) 
This is the mechanism for pressure gradient-trapping. 
Places with a negative value of V • F should produce 
a high particle density (see JK05). The correlation be- 
tween V • F and n is existent, but is very weak. The last 
panel, however, shows that there is a clear correlation 
between gas density and particle density, as is also evi- 
dent from Fig. 2] All in all, the correlations, even though 
some of the are quite weak, give the necessary informa- 
tion about the source of the dust concentrations. The 
concentrations are primarily due to pressure gradient- 
trapping in the gas flow. There is also evidence of some 
vorticity-trapping happening on top of that. 

Increases in density of up to two orders of magnitude 
will make a difference in the coagulation process, because 
at places of larger concentration more collisions (both 
destructive and constructive) are possible. Also there is 
a chance of increasing the density to such high values that 
a gravitational instability can occur in the densest places. 
We will consider this last point in more detail in Sect. [3] 
In the following section we show that the turbulence not 
only causes concentrations, but also changes the radial 
drift velocity of the boulders. 




Fig. 5. — Distribution of maximum particle densities. The curves 
show the fraction £(n) of particles that have been part of a given 
particle density during the 100 orbits. For i?o T f = 0.1, only con- 
centrations up to 10 are common, whereas for OqT{ = 1, 70% of 
the particles have experienced at least a 10 times increase in den- 
sity and 2% even a 20 times increase. For massive boulders with 
^0 T f = 10, more than 10% were part of a 30 times increase in 
density. 

ping of particles is happening in anticyclonic regions, 
and that regions of cycl onic flow are expelling particles 
ijBarge fc Sommerialll995lh Meter-sized particles should 
be optimally concentrated by vorticity, so the weak cor- 
relation between n and w z is surprising, considering that 
for centimeter-sized particles, JK05 find an almost linear 
relation between n and oj z with very small spread. The 
explanation may be that the friction time is so high that 
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Fig. 6. — Correlations between particle number density and various gas parameters. The first row considers two components of the 
vorticity. There is some correlation between n and uj z , indicating that particles are trapped in regions of anticyclonic flow. In the second 
row we consider the correlation between particle density and pressure gradient flux (explained in the text) and gas density, respectively. 
The first correlation is very weak, whereas there is evidently a correlation between gas density and particle density, although the fluctuation 
bars are significant. 

Here we have weighted the drift velocity with the num- 
ber density so that we are effectively measuring the aver- 
age momentum. We consider now for simplicity particles 
that have been accelerated by the gas to their terminal 
velocity (eq. |15j including the fluctuation pressure gra- 
dient), 

P H t^- ) , (17) 



The global pressure gradient on the gas forces solids to 
fall radially inwards. If the gas motion in the disc was 
completely non-turbulent, then the equilibrium radial 
drift velocity arising from the head wind term present 
in equation Ijl2(l would be 



6 



fl Tf + (QoTf)~ 



(15) 



dx 



We derived this expression by solving for dv/dt = in 
equation (|12|) . The highest drift speed occurs for par- 
ticles with ]?o T f = 1 with a laminar drift velocity of 
v x /c s = [3/2. We have checked by putting particles of 
different friction times into a non-turbulent disc that the 
measured drift velocities are in complete agreement with 
equation Q15[l. 

The effect of a real turbulent disc on the average drift 
velocity is seen in Fig. Here the average radial veloc- 
ity of all the particles is shown as a function of time for 
run A. For reference we overplot the laminar drift ve- 
locity (v x = —0.02 c s ) from equation l|15|l and the time- 
averaged drift velocity (Wx~ = — 0.012 c s ). The mean drift 
velocity is noticeably affected by the turbulence and its 
absolute value is reduced by 40% compared to the lami- 
nar value. The influence that turbulence can have on the 
mean drift velocity of the particles can be quantified with 
some simple analytical considerations. Considering the 
particles for a moment as a fluid with a number density 
scalar field n and a velocity vector field w, the average 
radial velocity can be calculated with the expression 



where e is defined as e = l/[S7oTf + (i?oTf) -1 ]. Inserting 
now equation (|17|) into equation , the resulting drift 
velocity is found to consist of two terms, 



= e/3c s 



^dx 

OX 



(n)L x 



(18) 



f Xl nw x dx 
(n)L x 



(16) 



The first term on the right-hand-side of equation JTJ 
represents the contribution to the average drift velocity 
from the global pressure gradient (eq.^fj])- The other 
term is an extra contribution due to any non-zero cor- 
relation between number density n and radial pressure 
gradient dlnp/dx. This situation is sketched in Fig. [SJ 
Here we sketch the global density gradient (3 (full line) 
and a sinusoidal density fluctuation In p{x) (dotted line). 
Particles concentrate in regions where the gas density 
fluctuation is positive, because there the divergence of 
the particle velocity is negative. Due to the total pres- 
sure gradient, the newly produced particle clumps drift 
inwards until the point where the outwards drift towards 
the fluctuation density maximum balances the inwards 
drift from the global pressure gradient. This is exactly 
around the location of the box in Fig. El Here the correla- 
tion between n and d In p/dx leads to a positive value of 
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Fig. 7. — Average radial particle velocity as a function of time for meter-sized boulders. The non-turbulent drift velocity is v& 
— 0.02c s (indicated with a dashed line), while the average drift velocity in the turbulent case is only around = — 0.012 c s , a reduction by 
around 40% in speed. 



the integral in equation i|18|). A closer inspection of Fig. 0] 
reveals that the dust overdensities are situated slightly 
downstream of the gas density fluctuation peaks, which 
is in good agreement with the the prediction in Fig. |S| 
If a significant fraction of the particles end up in such 
regions, the average drift speed is reduced 2 . For runs B 
and C, there is no significant reduction of the drift speed 
(see Table EJ), but there the predicted drift speed is also 
ten times lower than for meter-sized objects. Thus the 
measurement is not as reliable because the random veloc- 
ity fluctuations of the particles dominate over the radial 
drift. 

Due to the periodic boundary conditions in the 
y-directions, density structures quickly pass the y- 
boundaries, by shear advection, and thus possibly have 
some interference with themselves. To see the effect of 
the toroidal box size on the radial drift, we have run sim- 
ulations with a box size of 1.32 x 5.28 x 1.32 (run E) and 
1.32 x 10.56 x 1.32 (run F), keeping the resolution con- 
stant by adding the appropriate number of grid points 
in the y-direction. The time evolution of the mean ra- 
dial drift velocity is shown in Fig. |5J It is evidently very 
similar to Fig. [TJ so the toroidal size of the box does not 
influence the radial drift reduction noticeably. As seen in 
Table |21 the maximum particle density for runs E and F 
is quite high at 50 times the average density in the box, 
but not as high as in run A. However, simulations E and 
F only ran for 24 and 16 orbits, respectively, because of 
computational requirements due to the many grid points. 

In simulations of the interaction between a 
planet and a mag n etorot ationally turbulent disc, 
iNelson fc Papaloizoul l)2004|) find that the average 
migration velocity of the planets is not changed by the 
presence of MRI turbulence (whereas the spread in drift 
velocity causes some planets to even d rift outwards) . 
On the other hand, recent simulations bv iNelsonl l)2005fl 
indicate that the mean migration of planets can indeed 
change because of turbulence. The fluctuations in 

2 A more graphic explanation of the speed reduction is to con- 
sider a car race over a distance of 100 km. Half of the distance 
is sand, where the cars can run 50 kilometers per hour, and the 
other half asphalt, where the cars go 150 kilometers an hour. The 
average speed of a single car reaching the finish line is less than 
100 kilometers per hour, simply because that car spent more time 
on sandy terrain than on asphalt. 



migration speed are however much stronger than the 
average (so that hundreds of orbits are needed for 
a reliable estimate of the average). This is a very 
different kind of drift behavior than for the boulders 
in the current work, where the fluctuations in the drift 
speed are actually much smaller than the average. The 
presence of long-lived attracting regions in the gas may 
be the reason why boulders react on turbulence in a 
completely different way than planets do. 

Diminishing the radial drift for meter-sized objects by 
roughly one half may not be saving the boulders from 
their fate of decaying into the star. One will have to 
investigate this process by additionally looking at the 
growth behavior of the boulders which are sweeping up 
small grains on their way inwards. This sweeping up is 
determined by the actual drift speed with respect to the 
local gas motion. Even if the mean drift speed is above 
the threshold for effective sticking, there will be phases 
of much lower radial drift, where growth can occur. The 
overdense regions would also greatly increase the rate of 
destructive encounters between larger bodies, and thus 
the reservoir of small bodies would be stronger replen- 
ished there. This would not only influence growth of 
larger bodies, but also possibly have observational con- 
sequences. 

The present simulations are done in the gentle situa- 
tion of turbulence in a local box. Global disc simulations 
have stronger turbulence and larger density fluctuations. 
One can predict that it would thus also lead to a larger 
decrease in radial drift speed. This would possibly give 
the meter-sized boulders enough time to grow to a size 
safe for radial drift. However, this yet has to be demon- 
strated in global simulations 3 . 

5. GRAVITATIONAL INSTABILITY 

We already showed that turbulence can strongly influ- 
ence the growth of boulders by slowing them down and by 
concentrating them locally. These results can be incor- 
porated into standard evolution codes for the solid ma- 

3 We have recently become aware of work done by 
Fromane & Nelson (2005) where the dynamics of boulders in mag- 
netorotational turbulence is considered in global simulations of ac- 
cretion discs. They found indeed that solids can be trapped inside 
persistent flow features for even a hundred orbits, i.e. the entire 
simulation length. 
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terial (e.g. lWeidenschillindll99fllDullemond fe Dominild 
2005), which try to grow planctcsimals from dust grains 
via coagulation. On the other hand the high local con- 
centration can also lead to a different way of planetesi- 
mal formation, i.e. gelation. In the gelation case a cloud 
of boulders is so dense that gravitational attraction be- 
comes important. While we will not study self-gravity 
by an A^-body approach in this work (as one should), 
we want at least demonstrate by simple estimations un- 
der what conditions the concentration of boulders could 
clump into planetesimals. 

The gravity constant G enters in self-gravity calcula- 
tions, and thus the equations are no longer scale-free, 
but depend on the adopted disc model. We characterize 
a disc model by a column density So, an average dust-to- 
gas mass density ratio eo (for boulders of the considered 
size range) and a scale-height-to-radius ratio of £. Of 
course, eo will be smaller than the global dust-to-gas ra- 
tio ~ 0.02, because only a part of the mass will be present 
in boulders of the considered size range. We choose for 
simplicity the value eo = 0.01, assuming that 50% of the 
total dust mass is in bodies of the considered size, and 
we shall later discuss in how far this value is reasonable 
for a protoplanetary disc. 

The apparently large number of particles in our numer- 
ical simulations is still orders of magnitude away from 
any real number of boulders in the volume of the pro- 
toplanetary disc considered in our simulations. Thus it 
is necessary and validated to let one superparticle repre- 
sent an entire swarm of many particles of similar location 
and velocity in the disc. Superparticle means in this con- 
text that one particle has the aerodynamic behavior of a 
single boulder, but represents a mass of trillions of such 
bodies as it mimics an entire swarm of protoplanetesi- 
mals. Similar assumptions are common in simulations of 
giant planet core formation from colli ding planetesimals 
(|Kokubo fc Idal200arThommes et~al]l2003f) as well as in 
cosmological TV-body simulations l)Sommer-Larsen et alJ 
12003^ . We let the simulation box represent the proto- 
planetary disc in the mid-plane. Each superparticle then 
contains the mass m = e\piV/N, where V is the volume 
of the box, N is the number of superparticles, and ei 
and pi are the dust-to-gas ratio and the gas density in 
the mid-plane of the disc. We shall use the isothermal 
disc expression p\ = Sq/{^/2ttH) to calculate the mass 
density in the mid-plane. 

To calculate the dust-to-gas ratio in the mid-plane, ei, 
one needs to take into account the effect of vertical set- 
tling of solid material. Solids move in the direction of 
higher gas pressure. In the case of vertical stratification, 
that means that the boulders must sediment towards the 
mid-plane. An equilibrium is reached when the sedimen- 
tation is balanced by the turbulent diffusion, w ith dif- 
fusion coefficient D t l)Schrapler fc Hennind 1200 41. away 
from the mid-plane. This leads to a Gaus sian profile of 
the dust-to-gas ratio ijDubrulle et alJll995|) . 



e = eiexp[-z 2 /( 2 # e 2 



(19) 



with the dust-to-gas ratio scale height given by the ex- 
pression Hj = D t / (tiQq). The dust-to-gas ratio at z = 
is 
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Fig. 8. — Sketch of how turbulent density fluctuations can cause 
the average drift velocity to change. The full line shows the global 
density as a function of radial distance from the center of the box. 
On top of this we sketch a large-scale sinusoidal density fluctuation 
(dotted line) and the total density (dashed line). Dust particles are 
concentrated in the positive part of the fluctuation. At the same 
time the concentration drifts towards the location of the box where 
the total drift speed is zero. If a significant fraction of the particles 
end up in such regions, then the average drift speed can decrease. 
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Fig. 9. — Drift velocity for simulations E and F with larger y- 
domains. The expected drift velocity in a laminar disc is indicated 
with a dashed line. The measured drift velocity is approximately 
the same as for the cube simulations, so the periodic {/-boundary is 
not the reason for the reduced drift speed. Rather it is a side-effect 
of trapping the particles in radial density enhancements. 

where H = CgJ?^ 1 is the scale height of the gas. We 
now proceed by writing the turbulent diffusion coefficient 
as D t = <StCgi?Q~~ . where 5 t is the turbulen t diffusion 
equivalent of a t of lShakura fc Sunvaevl JT973). Then the 
mid-plane dust-to-gas ratio t\ can be written as 



i 



[2 T{ 



(21) 



(20) 



where the approximate expression is valid for i?o T f 3> <$t ■ 
For S t = a t = 0.002 and f2 Tf = 1, this gives ei ps 22.4e , 
so starting from a dust-to-gas ratio of eo = 0.01, the 
mid-plane dust-to-gas ratio can be expected to rise to 
ei = 0.22 due to vertical settling. Such a low dust- 
to-gas ratio alone will not for any physically reason- 
able column density cause gravitational fragmentation 



Gravoturbulent Formation of Planetesimals 



11 



fl Tpl.O , a t =0.002 , r =5.0AU , ff/r=0.04 , I =900.0 g cm" 2 , e =0.01 , e,=0.22 , £=50.9 orbits 




0.001 



Fig. 10. — Particle number density n in units of average density no, velocity dispersion <r in units of sound speed c s , free-fall time tff 
relative to the clump life-time t c \ , and clump radius R. together with Jeans radius Rj , all as a function of the number of included particles 
around the densest grid point in the box at a time of 50.9 orbits of run A. The vertical and horizontal dot-dot-dot-dashed lines indicate 
the regions of gravitational instability for the choice of disc model parameters. 



ijGoldreich fc Wardl 1 1973(1 or be subject to vertical stir- 
ring by the Kelvin-Hclmholtz instabilit y (the Rich ard- 
son number Ri is around unity, see e.g, ISekivalll998l and 
stratification with Ri > 0.25 should be stable). Even 
at such a high dust-to-gas ratio we are still in the gas- 
dominated regime where the back-reaction from the dust 
on the gas can be neglected. The turbulent dust concen- 
trations are assumed to occur in such a vertically settled 
dust layer. Now the most overdense regions will have 
a dust-to-gas ratio of unity and beyond. But we have 
measured that only about 3% of the grid cells have a 
dust-to-gas ratio of above unity at any given time, and 
thus it is still reasonable as a first approximation to ig- 
nore the back-reaction of the dust on the gas, although 
a more advanced study should include this effect as well. 

To find out if a given overdense clump is gravitationally 
unstable, we shall compare the different time-scales and 
length-scales involved in fragmentation by self-gravity in 
a Jeans- type stability analysis. First we investigate if the 
clump is gravitationally bound. We consider a clump 
of radius R, mass M and velocity dispersion a. The 
velocity dispersion must include the dispersion due to the 
background shear. For such a clump with a given mass 
to be gravitationally unstable, it must have a radius that 



is smaller than the Jeans radius given by 

2GM 

R j = 7T- ■ 



(22) 



If this first criterion is fulfilled, then it is also important 
that the collapse time-scale of the structure is shorter 
than the life-time of the overdense clump t c \. Only then 
we can be sure that the changing gas flow will not dissolve 
the concentration before it has had time to contract sig- 
nificantly. The fragmentational collapse happens on the 
free-fall time-scale 



tff = 



R 3 
GM 



(23) 



The condition for gravitational instability is now that 
R < Rj and that at the same time tg < t c \. We do 
not have to check separately that the collapse happens 
faster than a shear time t s h = S^q 1 , since the effect of 
the background shear is already included in the velocity 
dispersion. 

We now try to find out the smallest value of So that 
gives rise to a gravitational instability. Then we can see 
whether this is a value that occurs in nature or not. For 
J?oTf = 1 (run A) the minimum value of the column 
density turns out to be around Sq — 900 gem -2 (6 times 
the minimum mass solar nebula value at 5 AU) , whereas 



12 



Johansen et al. 



HgTplO.O , oc t =0.002 , r =5.0AU , #/r=0.04 , E =150.0 g cm" 2 , e =0.01 , £,=0.71 , £=53.0 orbits 
1000 1 \ i I i i 0.05 




0.001 



Fig. 11. — Same as Fig. E3 but for run C (decameter-sized bodies) at a time of 53 orbits. Here the minimum mass solar nebula column 
density is sufficient to have a gravitational instability. This is mainly because the velocity dispersion is smaller than for run A. Also the 
high density region has a larger extent. 



for f2oT{ = 10 (run C), a gravitationally unstable cluster 
of protoplanetesimals is achieved already at the minimum 
mass solar nebula value So = 150gcm~ 2 . 

First run A is considered. Here we can calculate the 
mass in each superparticle. With Sq = 900 gem -2 , 
£ = 0.04, r = 5 AU and ei = 0.22, we get Pl 
1 :2 x 10~ 10 gcm~ 3 and m = 8 x 10 20 g. Thus each su- 
perparticle represents about 3 x 10 14 meter-sized pro- 
toplanetesimals. This is five orders of magnitude more 
mass than in a kilometer-sized planetesimal, but since 
we are interested in identifying gravitationally bound re- 
gions with the mass of thousands and thousands of plan- 
etesimals, this is not a problem. Actually resolving the 
mass of even one single planetesimal with meter-sized 
objects would require on the order of a billion particles, 
which is way beyond current computational resources. 

We examine the region around the densest grid point 
of run A at a time of t = 50.9 orbits in more detail. This 
time is chosen because there occurs a large concentration 
of particles, see Fig. ^ We consider the j nearest par- 
ticles to the densest point and calculate for j between 1 
and 200,000 the particle number density n, the velocity 
dispersion er and its directional components, the free-fall 
time is , and the radius of the clump together with the 
Jeans radius Rj. The results are shown in Fig. It is 
reasonable to require at least j = 100 for a measurement 



to be statistically significant (for j > 100 the relative 
count ing error falls below 10%, see e.gT ICasertano fc Hut! 
1985). It is also reasonable to require that the size of 
the clump be larger than the size of a grid cell, since any 
structure in the concentration within a single grid is not 
well-resolved. The same is true for the velocity disper- 
sion. At j = 100 the dust number density is more than 
130 times the average, but the radius of the j = 100 
clump is only around 0.007, which is smaller than the 
grid cell radius of Sx/2 — 0.01. At j ~ 500 the clump 
has the size of a grid cell, and here the number density is 
more like 100 times the average. This must be multiplied 
by the enhancement by sedimentation, which is around 
20, to give a dust-to-gas ratio increase by a factor of 
2000 compared to the original value in the disc. The 
velocity dispersion is around a ~ 0.02 ... 0.03 c s . That 
includes the velocity dispersion due to the background 
shear, but this is not a very important effect anyway be- 
cause the size of the overdense clump is very small. At 
small scales the velocity dispersion is completely dom- 
inated by the radial component, according to Fig. ^3 
whereas the shear only takes over at larger scales. 

The free-fall time is a bit below the clump life-time, 
which is typically one shear time (see insert in Fig. ^ 
note that the time unit is in orbits). For calculating the 
Jeans radius we have had to adopt a column density as 
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Fig. 12. — Average velocity dispersion and fluctuation interval as a function of the number of particles in a grid cell. The dispersion 
rises until there are around 50 particles in a grid cell, and is then constant up to 200 particles, or around 30 times the average dust density. 
This corresponds to an isothermal equation of state for the boulders. 

high as Sq = 900 gem -2 in order to have the clump to be 
gravitationally unstable. This is mainly due to the high 
velocity dispersion. The radius of the clump is around 
one Jeans radius at j — 1000, so the clump is gravi- 
tationally bound at this scale and would be subject to 
further contraction by self-gravity. The gravitationally 
unstable region is around three grid cells in diameter, 
but even though this is well within the dissipative scales 
of the turbulence, the effect of the unresolved turbulence 
on the motion of the particles should be very little, as 
such small scale turbulence has short life-times and low 
amplitudes compared to the large scales. Extrapolating 
the resolved large scale turbulence to the grid-scale with 
a Kolmogorov law gives lower turbulent velocities than 
the particle velocity dispersion that we already measure 
at the grid scale. Thus we conclude that the unresolved 
turbulence has little or no influence on the particle dy- 
namics. The concentrations and velocity dispersions are 
exclusively driven by the large resolved scales of the gas 
motion. 

The solid size of the forming object would be roughly 
400 km if all the 1000 superparticles end up in just one 
large body. On the other hand, the outcome of such a 
collapse may also favor the further fragmentation of the 
clump. This all depends on how the velocity dispersion 
behave s with increasing de nsity. In the ./V-body simula- 
tions of lTanga et alJ l)2004j) gas drag works as an efficient 
way to dissipate the gravitational energy that is released 
in the contraction of protoplanetesimal clusters. Only 
such simulations, that include self-gravity and gas drag, 
could show the further evolution of the overdense boulder 
clumps that we see in the present work. 

For decameter-sized bodies (run C), we plot in Fig. ITT1 
the same quantities as in Fig. llOl around the densest point 
at a time of 53 orbits. This time we adopt the minimum 
mass solar nebula column density of Eq — 150 gem -2 , 
which gives a mid-plane density of pi = 2 x 10 -11 g cm" 3 . 



Because of the high friction time, the dust-to-gas ratio 
in the mid-plane fea.|21|) is now 0.71. The Richardson 
number is correspondingly lower at around Ri = 0.4, so it 
is still stable to Kelvin-Helmholtz instability. The mass 
of the individual superparticles is here m = 4 x 10 20 g. 
The density is slightly smaller than for meter-sized bod- 
ies, at statistically significant counts around 50 times the 
average, but the velocity dispersion is lower, and also the 
overdense region is much larger than it was for meter- 
sized boulders. Thus already a minimum mass solar neb- 
ula can produce a gravitational instability. The unstable 
region is as large as 10 grid cells in diameter, and contains 
around 10 5 particles. The size of a solid object consisting 
of this number of superparticles is roughly 1,400 kilome- 
ters. Again there is also the possibility that millions of 
10-kilomcter objects form instead. 

The preceding calculations are of course only an es- 
timation of the potential importance of self-gravity. In 
a real protoplanetary disc there will be a distribution 
of dust grain sizes present at any time. If e.g. frag- 
mentation is important, as discussed in the introduction, 
then the greater part (80%) of the mass may still be 
prese nt in bodies that are well be low one meter in ra- 
dius Ipullemond fc Dominikll2005D . With only 20% of 
the mass in the size range between one and ten meters, 
the critical column density could be as much as a factor 
of two higher than stated above. However, this is still in 
the range of the masses derived for circumstellar discs. 
So the qualitative picture that the clumps are gravita- 
tionally unstable for physically reasonable gas column 
densities is robust. 

To quantify the velocity dispersion in the entire box, 
we have calculated the average values over all the grid 
cells. The results are shown in the last four columns 
of Table [3 Grid cells with or 1 particles have been 
excluded from the average because the velocity disper- 
sion is per definition zero in these underresolved cells. 
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The meter-sized bodies have the highest velocity disper- 
sion, around o\ ~ 0.02 c s , whereas decimeter bodies have 
cq.i ~ 0.014 c s and decameter bodies have a value of 
trio w 0.017 c s . These values are similar to the turbulent 
velocities of the gas at the largest scales of the box (see 
Fig. 2 in JK05), which again shows that these large scales 
are the drivers of the particle dynamics. Interestingly run 
D, which is similar to run A only without the radial pres- 
sure gradient, has the same velocity dispersion as run A, 
so the radial pressure gradient does not add extra veloc- 
ity dispersion to the boulders. The toroidal component 
of the velocity dispersion is similar for all the runs be- 
cause it is dominated by the shear over a grid cell. Run 
C has a twice as large radial velocity dispersion as run 
B. This can be explained because the large particles in 
run C react much slower to the local behavior of the gas, 
and thus particles of different velocities and histories are 
mixed in together. 

The behavior of the velocity dispersion with increas- 
ing dust density is relevant for gravitational instability 
calculations. The average velocity dispersion, and the 
fluctuation width, as a function of the number of parti- 
cles in a grid cell is shown in Fig. El Again it is evident 
that the velocity dispersion for SloTf °f unity is largest. 
For all runs the velocity dispersion typically rises until 
there are around 50 particles in the cell. Then the dis- 
persions stay constant all the way to 200 particles. Thus 
the equation of state of the particles is isothermal, at 
least up to 30 times the average dust density. 

6. SUMMARY AND DISCUSSION 

We have considered the effect of magnetorotational 
turbulence on the motion of dust particles with a freely 
evolving space coordinate. The particle treatment was 
necessary over the fluid treatment, because the mean free 
path of the macroscopic dust boulders is so long that they 
can no longer be treated as a fluid. The use of magne- 
torotational turbulence may not be completely justified 
in the mid-plane of the disc where the ionization fraction 
due to radiation and cosmic particles is low. But due to 
its Kolmogorov-like properties, where energy is injected 
at the unstable large scales and then cascades down to 
smaller and smaller scales, magnetorotational turbulence 
can be seen as a sort of "generic disc turbulence" . 

We find that the turbulence acts on the particles by 
concentrating meter-sized boulders locally by up to a fac- 
tor of 100 and by reducing their radial drift by 40%. Both 
the concentrations and the reduced radial drift happen 
because the dust particles are temporarily trapped in ra- 
dial density enhancements. One would not expect such 
structures to be long-lived in a general turbulent flow, 
but magnetorotational turbulence in accretion discs is 
subject to a strong shear that favors elongated toroidal 
structures. In the presented simulations the typical life- 
time of the structures is on the order of a few orbits, 
corresponding to tens or even hundreds of years in the 
outer parts of a protoplanctary disc. When the density 
structures eventually dissolve, new structures appear at 
other locations. We find a strong correlation between a 
gas column density of a few percent above the average 
and a several times increase in the dust column density. 
We have also seen some evidence for increased dust den- 
sity in regions of anticyclonicity, but the long friction 
time of the dust particles makes it difficult to identify 



the gas flow that caused a given concentration, because 
the concentration may drift away from the creation site. 

The large concentrations naturally occur near the grid 
scale. In finite resolution computer simulations the dis- 
sipative length scale must necessarily be moved from the 
extremely small dissipative scales of nature to the small- 
est scales of the simulation box. Thus the turbulence 
is not well-resolved near the grid scale. On the other 
hand, the concentrations are driven by the largest scales 
of the turbulence, because there are the largest v eloci- 
ties and the longest lived features ijVolk et al.lll980() . Al- 
ready the other well-resolved but slightly smaller scales 
fluctuate too quick and at too low speeds to influence 
the path of an object that is one meter in size or larger. 
This argument is given support by the fact that we mea- 
sure particle velocity dispersions in the grid cells that are 
comparable to the velocity amplitude of the gas at the 
largest scales of the simulation. Thus, one should not 
expect higher resolution to change the concentrations or 
the velocity dispersions significantly. 

Our estimation of the minimum gas column density 
that would make the densest protoplanctesimal clumps 
gravitationally unstable is necessarily based on many as- 
sumptions. We assumed that half of the dust mass in the 
disc was present in bodies of the considered size, whereas 
in real discs an even larger part of the dust mass may be 
bound in small fragments that result from catastrophic 
collisions. We also ignored the back-reaction from the 
dust on the gas. The background state has, both for me- 
ter and decameter bodies, a dust-to-gas ratio just below 
unity (where the back-reaction becomes important). The 
effect of dust drag on the magnetorotational instability 
has to our knowledge never been considered. One can 
speculate that the drag force will mimic a strong vis- 
cosity and thus disable the source of turbulence where 
the dust density is high. For the treatment of Kelvin- 
Helmholtz instability we based it simply on a criterion 
on the Richardson number Ri. There is some indication 
that this may be too simplistic and that in protoplanc- 
tary di scs much higher Richards on numbers are also un- 
stable l|Q6mez fc Ostrikerll2005l) . but one can also spec- 
ulate that the full inclusion of dust particles in simula- 
tions of Kclvin-Helmholtz turbulence would show strong 
local concentrations like we see here for magnetorota- 
tional turbulence. Thus the exact values of six times the 
minimum solar nebula for meter-sized boulders and just 
the minimum mass solar nebula for decameter-sized boul- 
ders should only be considered as rough estimates. Still, 
the result that the clumps are gravitationally unstable 
for reasonable gas column densities is robust enough to 
warrant further investigations that include treatment of 
self-gravity between the boulders. 

Thus we find that the gravoturbulent formation of 
planetesimals from the fragmentation of an overdense 
swarm of meter-sized rocks is possible. Turbulence is 
in this picture not an obstacle, but rather the ignition 
spark, as it is responsible for generating the local gravita- 
tionally bound overdensities in the vertically sedimented 
layer of boulders. 



Computer simulations were performed at the Danish 
Center for Scientific Computing in Odense and at the 
RIO cluster at the Rechenzentrum Garching. Our re- 
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